The Dataset

The observations in the dataset milan.mort were made to check whether the presence of small particles (namely SO2 and TSP) can be associated with a higher number of deaths. The study was conducted in the city of Milan, by Vigotti, M.A., Rossi, G., Bisanti, L., Zanobetti, A. and Schwartz, J. in a paper called “Short term effect of urban air pollution on respiratory health in Milan, Italy, 1980-1989”.

The dataset has 9 variables:

Each of these quantities has been observed every day for ten years (from 1980 to 1989). This amounts to 3652 observations.

Our goal is to model the time series for the dependent variable tot.mort using both a Gaussian and a Poissoinian likelihood. We will try to use every explanatory variables except resp.mort (since the information carried by this variable is somehow already included in tot.mort).

Exploratory Data Analysis

First off, we add two more varibales:

data$day.of.year = data$day.num %% 365
data[data$day.of.year == 0, 'day.of.year'] = 365
data$year = data$day.num %/% 365

Next, we plot the evolution through time of both the total and respiratory number of deaths.

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

The data are extremelly noisy, but thank to the local regression (loess) smoothing we can spot a cyclic pattern which is repeted across the ten years.

We now plot the distribution of the dependent variables:

We can see that SO2 is strongly skewed suggesting that logarithmic transformation is necessary:

data$SO2 <- log(data$SO2 + 1 + abs(min(data$SO2))) #traslation for negative values

On the other hand mean.temp shows a quite uniform distribution and tot.mort follows a centered bell-shaped distribution. Also we can see that the number of resp.mort is very low compared to tot.mort.

Finally, we can look at the pairplot for more insights:

cor.data <- cor.matrix(variables=d(mean.temp,rel.humid,SO2, TSP,tot.mort),
                                   data = data)

ggcorplot(
  cor.mat = cor.data,
  data = data,
  line.method=c("loess"),
  var_text_size = 3,
  cor_text_limits = c(5,10))
## `geom_smooth()` using formula 'y ~ x'

We can see that the effect of the independent variables is kind of weak. This suggest narrow prior centered at zero for the slopes. Also, the smoothing highlights some nonlinearity, especially in mean.temp, for which we can see a minimum in the number of deaths for mild temperatures.

One last assumption we are making is that holiday and day.of.week do not affect the final result. Indeed we assume an homogeneuos distribution for tot.mort over day.of.week and the same ratio of deaths on any working day or holiday. Indeed we have

ggplot(data=data, aes(x = as.factor(day.of.week), y=tot.mort, fill=as.factor(day.of.week), color = as.factor(day.of.week)))+geom_bar(stat = "identity")+ggtitle(label = "Weekday Mortality Comparison")

and

n_holidays <- sum(data[data$holiday==1,'holiday'])
print(sum(data[data$holiday==1, 'tot.mort'])/n_holidays)
## [1] 31.71569
print(sum(data[data$holiday==0, 'tot.mort'])/(nrow(data)- n_holidays))
## [1] 31.85904

Note that we exclude also a weekend effect because we are considering a small sample (only the city of Milan) and the impact factors such of work accidents or traffic accident is negligible.

Gaussian Models

We now build a some models supposing a normally distributed response variable.

Initial Simple Linear Model (‘G0’)

We start with a linear model which does not explicitly include the effect time (implicitly this information may be contained in mean.temp). We will use this model as the baseline.

As mentioned before, we will set a gaussian with small variance for the parameters; on the other hand we will give the intercept more freedom as it may be harder to guess. Therefore the model will be

\[ \text{tot.mort} \sim \mathcal{N}(\alpha + \beta_1 \text{temp} + \beta_2 \log\text{SO2}, \sigma^2) \]

with the following priors

\[ \beta_i \sim \mathcal{N}(0,15)\\ \]

we will use the default weakly informative priors provided by rstanarm

priors <- normal(location = c(0,0), scale=c(15,15), autoscale = FALSE) #prior for the slopes

fitted.gaussian_regression <- stan_glm(tot.mort ~ mean.temp + SO2, data = data, family = gaussian, prior=priors )
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 8.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.85 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.210172 seconds (Warm-up)
## Chain 1:                0.493138 seconds (Sampling)
## Chain 1:                0.70331 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 2.5e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.265892 seconds (Warm-up)
## Chain 2:                0.464404 seconds (Sampling)
## Chain 2:                0.730296 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.209692 seconds (Warm-up)
## Chain 3:                0.47755 seconds (Sampling)
## Chain 3:                0.687242 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 2.2e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.199602 seconds (Warm-up)
## Chain 4:                0.525765 seconds (Sampling)
## Chain 4:                0.725367 seconds (Total)
## Chain 4:

We now make use replicated data to perform some checks.

y_rep <- posterior_predict(fitted.gaussian_regression)

First of all we check wheter the residuals are normally distributed

mean_y_rep <- colMeans(y_rep)
std_resid <- (data$tot.mort - mean_y_rep) / sqrt(mean_y_rep)

qplot(mean_y_rep, std_resid) + hline_at(2) + hline_at(-2)

Some large residuals are located at mean_y_rep\(\simeq 27\), but overall they seem to be randomly distributed.

Next we check on residual autocorrelation

acf(data$tot.mort - mean_y_rep)

We see that The model may not be correctly specified as there is large poitive autocorellation at all lags.

We now proceed with the other usual posterior predictive checks:

ppc_dens_overlay(
  y = sapply(data$tot.mort,as.numeric),
  yrep = y_rep[1:200,]
)

ppc_stat_grouped(
  y = sapply(data$tot.mort,as.numeric), 
  yrep = y_rep, 
  group = data$year, 
  stat = 'mean',
  binwidth = 0.2
)

We can see that the model does not perform very well: it doesn’t reproduce the empirical distribution as the peak is somewhat shifted and the mean value is completely wrong in different years.

The LOOIC is

loo_slopes1 <- loo(fitted.gaussian_regression)
print(loo_slopes1)
## 
## Computed from 4000 by 3649 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -12129.4  53.9
## p_loo         4.8   0.4
## looic     24258.9 107.8
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

We have seen the residuals and the mean in the ppc_stat_grouped plot are the two main issues with the model. Therefore we will use these together with the LOOIC as benchmarks through which test the performance of the following models.

Gaussian Regression with Time

As it suggested from EDA, the time period can be important in predicting the total number of deaths. So we now explore some models which include time effect.

First Model Including Time (‘G1’)

The model specification are the following:

\[ \mu = \alpha + \beta_1 \text{day.of.year} + \beta_2\text{year} + \beta_3 \text{temp} \\ \text{tot.mort} \sim \mathcal{N}(\mu, \sigma^2) \]

For this model we will use the default weakly informative priors of rstanarm.

fitted.time_regression <- stan_glm(tot.mort ~ year + day.of.year + mean.temp, data = data, family = gaussian)
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.7e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.090515 seconds (Warm-up)
## Chain 1:                0.414339 seconds (Sampling)
## Chain 1:                0.504854 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 2.2e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.115938 seconds (Warm-up)
## Chain 2:                0.413201 seconds (Sampling)
## Chain 2:                0.529139 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 2.1e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.21 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.093901 seconds (Warm-up)
## Chain 3:                0.409352 seconds (Sampling)
## Chain 3:                0.503253 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 2.1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.21 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.091321 seconds (Warm-up)
## Chain 4:                0.414289 seconds (Sampling)
## Chain 4:                0.50561 seconds (Total)
## Chain 4:

Again, we first look at the residuals

y_rep2 <- posterior_predict(fitted.time_regression)
mean_y_rep <- colMeans(y_rep2)
std_resid <- (data$tot.mort - mean_y_rep) / sqrt(mean_y_rep)
qplot(mean_y_rep, std_resid) + hline_at(2) + hline_at(-2)

acf(data$tot.mort - mean_y_rep)

Where we see that the residuals still suffer from positive autocorrelation but to a lesser degree with respect to the base model.

And then the posterior predictive checks:

ppc_dens_overlay(
  y = sapply(data$tot.mort, as.numeric),
  yrep = y_rep2[1:200,]
)

We can see that the replicated distribution somewhat replicate the original one, but the peaks are not aligned.

ppc_stat_grouped(
  y = sapply(data$tot.mort, as.numeric), 
  yrep = y_rep2, 
  group = as.factor(data$year), 
  stat = 'mean',
  binwidth = 0.2
)

We can see how modelling time has helped improving both the autocorrelation and the mean estimate for each year. However we are still not satisfied with this model

Finally, the LOOIC:

loo_slopes2 <- loo(fitted.time_regression)
loo_slopes2
## 
## Computed from 4000 by 3649 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -12073.1  53.7
## p_loo         5.7   0.4
## looic     24146.2 107.4
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Best Gaussian Model with Time (‘G2’)

After some tests, we saw that only one predictor leads to a significant reduction of LOOIC: log(SO2). Adding more predictors does not carry any significant benefits and we would simply end up in having a more complex model.

fitted.time_s02_regression <- stan_glm(tot.mort~year + day.of.year + mean.temp + SO2, data = data, family = gaussian)
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.135085 seconds (Warm-up)
## Chain 1:                0.468478 seconds (Sampling)
## Chain 1:                0.603563 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 2.1e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.21 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.134481 seconds (Warm-up)
## Chain 2:                0.486326 seconds (Sampling)
## Chain 2:                0.620807 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 2.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.182078 seconds (Warm-up)
## Chain 3:                0.463384 seconds (Sampling)
## Chain 3:                0.645462 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 2.1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.21 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.185112 seconds (Warm-up)
## Chain 4:                0.4685 seconds (Sampling)
## Chain 4:                0.653612 seconds (Total)
## Chain 4:

We perform the same steps as before:

y_rep3 <- posterior_predict(fitted.time_s02_regression)
mean_y_rep <- colMeans(y_rep3)
std_resid <- (data$tot.mort - mean_y_rep) / sqrt(mean_y_rep)
qplot(mean_y_rep, std_resid) + hline_at(2) + hline_at(-2)

acf(data$tot.mort - mean_y_rep)

#Result evaluation
ppc_dens_overlay(
  y = sapply(data$tot.mort, as.numeric),
  yrep = y_rep3[1:200,]
)

ppc_stat_grouped(
  y =  sapply(data$tot.mort, as.numeric), 
  yrep = y_rep3, 
  group = as.factor(data$year), 
  stat = 'mean',
  binwidth = 0.2
)

loo_slopes3 <- loo(fitted.time_s02_regression)
loo_slopes3
## 
## Computed from 4000 by 3649 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -12022.4  54.0
## p_loo         6.7   0.4
## looic     24044.7 107.9
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

As we can see, the only improvement we achieved is a reduction in LOOIC. Indeed, the posterior predictive check do not highlight any other significant achievement.

Nonlinear Gaussian Model (‘G3’)

EDA suggested that some predictors may be nonlinear in predicting the total number of deaths. We therefore build a gam model to try to take into account these effects. After some attempts, we found that the minimal best model is the following

library(rstanarm)

fitted.nonlinear <- stan_gamm4(tot.mort ~ s(mean.temp) + s(day.of.year) + factor(year) + SO2,
                                 family = gaussian, data = data)
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 4.7e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.47 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 8.30233 seconds (Warm-up)
## Chain 1:                6.02256 seconds (Sampling)
## Chain 1:                14.3249 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 3.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.34 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 7.79747 seconds (Warm-up)
## Chain 2:                9.04056 seconds (Sampling)
## Chain 2:                16.838 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 3.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.32 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 8.33332 seconds (Warm-up)
## Chain 3:                7.22921 seconds (Sampling)
## Chain 3:                15.5625 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 3.4e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.34 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 8.02176 seconds (Warm-up)
## Chain 4:                9.48606 seconds (Sampling)
## Chain 4:                17.5078 seconds (Total)
## Chain 4:
y_rep4 <- posterior_predict(fitted.nonlinear)
mean_y_rep <- colMeans(y_rep4)
std_resid <- (data$tot.mort - mean_y_rep) / sqrt(mean_y_rep)
qplot(mean_y_rep, std_resid) + hline_at(2) + hline_at(-2)

acf(data$tot.mort - mean_y_rep)

We see that the autocorrelation is much smaller, but it still positive up to the 15th lag.

ppc_dens_overlay(
  y = sapply(data$tot.mort, as.numeric),
  yrep = y_rep4[1:200,]
)

Again, the peaks are not aligned.

ppc_stat_grouped(
  y = sapply(data$tot.mort, as.numeric), 
  yrep = y_rep4, 
  group = as.factor(data$year), 
  stat = 'mean',
  binwidth = 0.2
)

This model finally gets the distribution of the mean centered. Let’s see how it behaves with the standard deviation distributions:

ppc_stat_grouped(
  y = sapply(data$tot.mort, as.numeric), 
  yrep = y_rep4, 
  group = as.factor(data$year), 
  stat = 'sd',
  binwidth = 0.2
)

loo_slopes4 <- loo(fitted.nonlinear)
loo_slopes4
## 
## Computed from 4000 by 3649 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo -11788.4 47.1
## p_loo        31.8  1.6
## looic     23576.9 94.1
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Up to now this is the best model: distribution of the mean is nicely cetered for each year and it achieves the smallest LOOIC. However, the replicated distributions and the residuals still show some issues.

Poissonian Models

Poisson Regression with Time

We now build some models with a Poisson distribution as likelihood. We expect this to perform better (in terms of LOOIC) as it should better represent count data. However we also expect that it will not solve the main issue of positive autocorrelation in the residuals.

First Model including Time (‘P1’)

Again, we use the default weakly informative priors of rstanarm.

fitted.P_time_regression <- stan_glm(tot.mort ~ year + day.of.year + mean.temp, data = data, family = poisson)
## 
## SAMPLING FOR MODEL 'count' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000225 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.25 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 1.71325 seconds (Warm-up)
## Chain 1:                1.93895 seconds (Sampling)
## Chain 1:                3.65219 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'count' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000198 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.98 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.82189 seconds (Warm-up)
## Chain 2:                1.99666 seconds (Sampling)
## Chain 2:                3.81855 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'count' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000198 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.98 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 1.70195 seconds (Warm-up)
## Chain 3:                2.06358 seconds (Sampling)
## Chain 3:                3.76553 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'count' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000194 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.94 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 1.72567 seconds (Warm-up)
## Chain 4:                2.08691 seconds (Sampling)
## Chain 4:                3.81257 seconds (Total)
## Chain 4:

Now the residuals

y_rep2P <- posterior_predict(fitted.P_time_regression)
mean_y_rep <- colMeans(y_rep2P)
std_resid <- (data$tot.mort - mean_y_rep) / sqrt(mean_y_rep)
qplot(mean_y_rep, std_resid) + hline_at(2) + hline_at(-2)

acf(data$tot.mort - mean_y_rep)

We again find positive autocorrelation in accordance with our expectation

#Result evaluation
ppc_dens_overlay(
  y = sapply(data$tot.mort, as.numeric),
  yrep = y_rep2P[1:200,]
)

We can see that the posterior predictive distributions better replicates the empirical one.

ppc_stat_grouped(
  y = sapply(data$tot.mort, as.numeric), 
  yrep = y_rep2P, 
  group = as.factor(data$year), 
  stat = 'mean',
  binwidth = 0.2
)

Here we can see that the distribution is totally wrong for year 2,3,4 and 6.

The model score is:

loo_slopes2P <- loo(fitted.P_time_regression)
loo_slopes2P
## 
## Computed from 4000 by 3649 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -12137.1  69.2
## p_loo         5.6   0.2
## looic     24274.2 138.5
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Expect for the ppc_dens_overlay plot, the model is worse than it gaussian counterpart

Best Poisson Model with Time (‘P2’)

Also in the poissonian models, the log(SO2) predictos turns out to be the best in terms of LOOIC improvement.

fitted.P_time_s02_regression <- stan_glm(tot.mort~year + day.of.year + mean.temp + SO2, data = data, family = poisson)
## 
## SAMPLING FOR MODEL 'count' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000198 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.98 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 2.36586 seconds (Warm-up)
## Chain 1:                2.513 seconds (Sampling)
## Chain 1:                4.87886 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'count' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000192 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.92 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 2.40657 seconds (Warm-up)
## Chain 2:                2.69719 seconds (Sampling)
## Chain 2:                5.10377 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'count' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000192 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.92 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 2.59173 seconds (Warm-up)
## Chain 3:                2.36845 seconds (Sampling)
## Chain 3:                4.96019 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'count' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000226 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2.26 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 2.22154 seconds (Warm-up)
## Chain 4:                2.5354 seconds (Sampling)
## Chain 4:                4.75694 seconds (Total)
## Chain 4:

And now the usual model checking steps:

y_rep3P <- posterior_predict(fitted.P_time_s02_regression)
mean_y_rep <- colMeans(y_rep3P)
std_resid <- (data$tot.mort - mean_y_rep) / sqrt(mean_y_rep)
qplot(mean_y_rep, std_resid) + hline_at(2) + hline_at(-2)

acf(data$tot.mort - mean_y_rep)

#Result evaluation
ppc_dens_overlay(
  y = sapply(data$tot.mort, as.numeric),
  yrep = y_rep3P[1:200,]
)

ppc_stat_grouped(
  y =  sapply(data$tot.mort, as.numeric), 
  yrep = y_rep3P, 
  group = as.factor(data$year), 
  stat = 'mean',
  binwidth = 0.2
)

loo_slopes3 <- loo(fitted.P_time_s02_regression)
loo_slopes3
## 
## Computed from 4000 by 3649 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -12078.9  67.8
## p_loo         6.8   0.2
## looic     24157.8 135.6
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Attempts for Improvement

In contradiction with our expectations, the poissonian models gave worse results with respect to the gaussian ones. This is probably due to the fact that a Poisson model restrics it variance to be equal to the mean. Indeed, such assumption may not be suitable for these noisy data.

We will therefore focus on gaussian models and try to improve them.

Hierarchial Gaussian Model (‘GH’)

As a first attempt we build a hierarchial model where year constitute the deeper level. This should allow us to somewhat model the yearly pattern of the number of deaths. This should also take into account the slightly deacresing trend in the number of deaths:

## `geom_smooth()` using formula 'y ~ x'

We also note that the variable SO2 assumes different mean values every year:

## `geom_smooth()` using formula 'y ~ x'

and therefore we will use a varying slope for this variable.

In conclusion, the model specifications are:

\[ \text{tot.mort}_{yi} \sim \mathcal{N}(\mu_{y(i)} + \beta_1 \text{day}_i + \beta_2 \text{temp}_i + \kappa_{y(i)}\text{SO2}_i, \sigma)\\ \mu_{y(i)} \sim \mathcal{N}(\alpha, \sigma_{\mu})\\ \kappa_{y(i)} \sim \mathcal{N}(\beta_{SO2}, \sigma_{\kappa}) \]

Note that we used a transforemed parameters block in order to avoid divergent series.

gaussian_hier_regression <- stan_model("01GaussianHier.stan")
stan_data_hier <- list(
  N = nrow(data), 
  J = length(unique(data$year)),
  day= data$day.of.year,
  temp = data$mean.temp,
  tot_mort = data$tot.mort,
  so2 = data$SO2,
  year_idx = data$year+1
)

fitted.gaussian_hier_regression <- sampling(gaussian_hier_regression, data = stan_data_hier)
## 
## SAMPLING FOR MODEL '01GaussianHier' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.001134 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.34 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 87.5474 seconds (Warm-up)
## Chain 1:                60.4431 seconds (Sampling)
## Chain 1:                147.991 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '01GaussianHier' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000552 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 5.52 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 84.2104 seconds (Warm-up)
## Chain 2:                61.9131 seconds (Sampling)
## Chain 2:                146.123 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '01GaussianHier' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000549 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 5.49 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 87.5713 seconds (Warm-up)
## Chain 3:                69.4471 seconds (Sampling)
## Chain 3:                157.018 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '01GaussianHier' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000551 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 5.51 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 91.3572 seconds (Warm-up)
## Chain 4:                51.0008 seconds (Sampling)
## Chain 4:                142.358 seconds (Total)
## Chain 4:
## Warning: There were 5 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
Computational Checks
grid.arrange(a1,a2,a3,a4,a5,a6, nrow=3, ncol = 2)

Graphically can see that overall the chains have mixed quite well

rhats_att_def <- rhat(fitted.gaussian_hier_regression, pars = c("alpha", "beta_day", 'beta_so2', 'sigma_mu', 'sigma_kappa', 'mu_raw', 'kappa_raw'))
mcmc_rhat(rhats_att_def)

Also the \(\hat{R}\) mostly close to one, showing again that the chains have mixed well.

n_eff_ratios <- neff_ratio(fitted.gaussian_hier_regression, pars = c("alpha", "beta_day", 'beta_so2', 'sigma_mu', 'sigma_kappa', 'mu_raw', 'kappa_raw'))

mcmc_neff(n_eff_ratios, size = 2)

Result evaluation
y_rep_hier <- as.matrix(fitted.gaussian_hier_regression, pars = "y_rep")
mean_y_rep <- colMeans(y_rep_hier)

acf(data$tot.mort - mean_y_rep)

#Result evaluation
ppc_dens_overlay(
  y = stan_data_hier$tot_mort,
  yrep = y_rep_hier[1:200,]
)

ppc_stat_grouped(
  y = stan_data_hier$tot_mort, 
  yrep = y_rep_hier, 
  group = as.factor(data$year), 
  stat = 'mean',
  binwidth = 0.2
)

log_lik_slopes <- extract_log_lik(fitted.gaussian_hier_regression)
loo_slopes_hier <- loo(log_lik_slopes, 
                   r_eff=relative_eff(exp(log_lik_slopes),rep(c(1,2,3,4),each=1000)))
loo_slopes_hier
## 
## Computed from 4000 by 3649 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -11989.6  52.3
## p_loo        18.6   0.7
## looic     23979.2 104.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Now the distributions of ppc_stat_grouped are centered, but overall the model does not seem to improve that much.

We can see that the variances of the varying slope and intercept are significantly different from zero:

p1 <- mcmc_hist(
  as.matrix(fitted.gaussian_hier_regression, pars = "sigma_mu"),
  binwidth = 0.05
)

p2 <- mcmc_hist(
  as.matrix(fitted.gaussian_hier_regression, pars = "alpha"),
  binwidth = 0.05
)

grid.arrange(p1,p2, nrow=1)

p1 <- mcmc_hist(
  as.matrix(fitted.gaussian_hier_regression, pars = "sigma_kappa"),
  binwidth = 0.05
)

p2 <- mcmc_hist(
  as.matrix(fitted.gaussian_hier_regression, pars = "beta_so2"),
  binwidth = 0.05
)

grid.arrange(p1,p2, nrow=1)

This suggests us that the random slope slopes and interceptes are important to explain the variability in the model and therefore we try to keep improving this model.

Hierarchical Gaussian Model with Autoregressive Effect (‘GHAR’)

Up to now we have included time in our models in a naive way. We now try to include an autoregressive daily effect.

\[ \text{tot.mort}_{yt} \sim \mathcal{N}(\mu_{y} + \text{day}_t + \beta_2 \text{temp} + \kappa_{y}\text{SO2}, \sigma)\\ \text{day}_t = \rho \text{day}_{t-1} + \epsilon_t, \text{ } \epsilon_t \sim \mathcal{N(0,\sigma_{day})}\\ \rho \in [-1,1] \] Note that while in the previous model \(\text{day}\) was equal to day.of.year, now it is equal to day.num.

We also have

\[ \text{day}_1 \sim \mathcal{N}({0, {\sigma_{day}}}/{\sqrt{1 - \rho^2}}) \] and

\[ \rho_{raw} \in [0,1]\\ \rho = 2 \times \rho_{raw} -1\\ \rho_{raw} \sim \mathcal{Beta}(3,2) \]

gaussian_autor_regression <- stan_model("02GaussianHierAuto.stan")
stan_data_auto <- list(
  N = nrow(data), 
  J = length(unique(data$year)),
  D = max(data$day.num),
  day_idx = data$day.num,
  temp = data$mean.temp,
  tot_mort =  data$tot.mort,
  so2 = data$SO2,
  year_idx = data$year+1
)

fitted.gaussian_autor_regression <- sampling(gaussian_autor_regression, data = stan_data_auto)
## 
## SAMPLING FOR MODEL '02GaussianHierAuto' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.005709 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 57.09 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 191.548 seconds (Warm-up)
## Chain 1:                132.219 seconds (Sampling)
## Chain 1:                323.768 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '02GaussianHierAuto' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000803 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 8.03 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 199.948 seconds (Warm-up)
## Chain 2:                133.212 seconds (Sampling)
## Chain 2:                333.16 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '02GaussianHierAuto' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.001008 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 10.08 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 189.341 seconds (Warm-up)
## Chain 3:                129.519 seconds (Sampling)
## Chain 3:                318.861 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '02GaussianHierAuto' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.001369 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 13.69 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 194.532 seconds (Warm-up)
## Chain 4:                130.561 seconds (Sampling)
## Chain 4:                325.094 seconds (Total)
## Chain 4:
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
Computational Checks
grid.arrange(a1,a2,a3,a4,a5,a6, nrow=3, ncol = 2)

rhats_att_def <- rhat(fitted.gaussian_hier_regression, pars = c("alpha", "beta_day", 'beta_so2', 'sigma_mu', 'sigma_kappa', 'mu_raw', 'kappa_raw'))
mcmc_rhat(rhats_att_def)

n_eff_ratios <- neff_ratio(fitted.gaussian_hier_regression, pars = c("alpha", "beta_day", 'beta_so2', 'sigma_mu', 'sigma_kappa', 'mu_raw', 'kappa_raw'))

mcmc_neff(n_eff_ratios, size = 2)

Also here we do not have any convergence or efficiency issue.

Result evaluation
y_rep_autor <- as.matrix(fitted.gaussian_autor_regression, pars = "y_rep")
mean_y_rep <- colMeans(y_rep_autor)
acf(data$tot.mort - mean_y_rep)

ppc_dens_overlay(
  y = sapply(data$tot.mort, as.numeric),
  yrep = y_rep_autor[1:200,]
)

Here we check the mean

ppc_stat_grouped(
  y = sapply(data$tot.mort, as.numeric), 
  yrep = y_rep_autor, 
  group = as.factor(data$year), 
  stat = 'mean',
  binwidth = 0.2
)

and here the sd:

ppc_stat_grouped(
  y = sapply(data$tot.mort, as.numeric), 
  yrep = y_rep_autor, 
  group = as.factor(data$year), 
  stat = 'sd',
  binwidth = 0.2
)

these plots are better with respect to the ones of the GAM model.

log_lik_slopes <- extract_log_lik(fitted.gaussian_autor_regression)
loo_slopes_autor <- loo(log_lik_slopes, 
                   r_eff=relative_eff(exp(log_lik_slopes),rep(c(1,2,3,4),each=1000)))
## Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
loo_slopes_autor
## 
## Computed from 4000 by 3649 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo -11687.0 46.5
## p_loo       315.6  7.7
## looic     23374.1 93.0
## ------
## Monte Carlo SE of elpd_loo is 0.3.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     3645  99.9%   532       
##  (0.5, 0.7]   (ok)          4   0.1%   532       
##    (0.7, 1]   (bad)         0   0.0%   <NA>      
##    (1, Inf)   (very bad)    0   0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.

We can see that this model outperforms all the previus ones. It also manges to achieve centered distriution for both the standard deviation and the mean for every year. However, the pakes of the empirical distribution and that of the replicated distributiions are still somewhat misaligned.

Conclusions

The data are very noisy and quite hard to model. Throughout this analysis we however managed to find out some aspects about the generative process. First of all we saw how a Gaussian likelihood is much more effective in describing the outcome, in spite of the fact that a Poisson should more suitable when dealing with count variables. Second, we saw that mean.temp and SO2 have the most significant effect on the overall mortality, while TSP and rel.humid don’t. We also saw that, because of the ciclical intrinsic nature of the data, time is an important predictor.

In addition, the GAM model that mean.temp has a nonlinear effect which can be effectively described through the spline based smooth s().

Finally, we managed to solve the autocorrelation in the residuals including an autogression on each day.

The obvious next step is to include the smoothing function in the hierarchical autoregressive model. However it would be interesting to explore different generative processes, such as the Negative Binomial, and models that are by nature more suitable for time series, such as ARIMA. In any case, adding more predictors (age, number of smokers) should definitly help in understanding the mortality in the city of Milan.